Load packages
load data
# load UMI
umi.table <- read.table("/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895282_Hashtag-RNA.umi.txt.gz",
header = TRUE, row.names = 1)
umi.table <- as.sparse(umi.table)
umi.table <- CollapseSpeciesExpressionMatrix(object = umi.table)
# load HTO
hto.table <- read.csv(file = "/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895283_Hashtag-HTO-count.csv",
header = TRUE, row.names = 1)
remove.index <- which(rownames(hto.table) %in% c("bad_struct", "no_match", "total_reads"))
hto.table <- hto.table[-remove.index, ]
# load ADT
adt.table <- read.csv(file = "/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895284_Hashtag-ADT2-count.csv",
header = TRUE, row.names = 1)
remove.index <- which(rownames(adt.table) %in% c("bad_struct", "no_match", "total_reads"))
adt.table <- adt.table[-remove.index, ]# only keep cells with RNA, ADT and HTO
joint.bcs <- intersect(colnames(umi.table), colnames(hto.table))
joint.bcs <- intersect(colnames(adt.table), joint.bcs)
umi.table <- umi.table[, joint.bcs]
hto.table <- as.matrix(hto.table[, joint.bcs])
adt.table <- as.matrix(adt.table[, joint.bcs])
rownames(x = hto.table) <- gsub(pattern = "-.+", replacement = "", x = rownames(x = hto.table))
rownames(x = adt.table) <- gsub(pattern = "-.+", replacement = "", x = rownames(x = adt.table))Create Seurat object from RNA data
Add HTO and ADT data
Demultiplex cells based on HTO enrichment
pbmc <- dataNormalization(object = pbmc)
pbmc <- HTODemux(pbmc, assay = "HTO", positive.quantile = 0.99)
table(pbmc$HTO_classification.global)
Idents(pbmc) <- "HTO_maxID"
VlnPlot(pbmc, assay = "HTO", features = rownames(pbmc[["HTO"]]), ncol = 3, pt.size = 0.01)# t-SNE visualization
pbmc <- subset(pbmc, idents = "Negative", invert = TRUE)
# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = pbmc, assay = "HTO"))))
# Calculate tSNE embeddings with a distance matrix
pbmc <- RunTSNE(pbmc, distance.matrix = hto.dist.mtx, perplexity = 100)
DimPlot(pbmc)joint analysis using RNA expression and ADT singals
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "mean.var.plot") # directly use Seurat Function
pbmc <- dataScaling(object = pbmc)
# PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), verbose = FALSE) # directly use Seurat Function
pbmc <- jointDistance(object = pbmc)
pbmc <- tsneFromDistane(object = pbmc, assay = "All")
pbmc <- clusteringFromDistance(object = pbmc, assay = "All")
# save(pbmc,file = 'HTO.rData')visualization ADT vs RNA vs Joint
plots <- generateGridDimPlot(pbmc, legend = FALSE, reduction.prefix = "tsne_",
darkTheme = FALSE)
listPlot(object = plots)Heat map that show both ADT and RNA
HTO on joint map
phto1 <- DimPlot(pbmc, group.by = "HTO_classification.global", reduction = "tsne_joint") +
DarkTheme()
phto2 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_joint") + DarkTheme()
phto3 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_rna") + DarkTheme()
phto4 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_adt") + DarkTheme()
plot_grid(phto1, phto2, phto3, phto4, ncol = 2)phto1 <- DimPlot(pbmc, group.by = "HTO_classification.global", reduction = "tsne_rna") +
xlab("t-SNE1") + ylab("t-SNE2")
phto2 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_rna") + xlab("t-SNE1") +
ylab("t-SNE2")
plot_grid(phto1, phto2, ncol = 2)features.adt <- rownames(pbmc@assays[["ADT"]])
features.hto <- rownames(pbmc@assays[["HTO"]])
FeaturePlot(pbmc, features = features.adt, min.cutoff = "q05", max.cutoff = "q95",
ncol = 3, reduction = "tsne_adt")FeaturePlot(pbmc, features = features.hto, min.cutoff = "q05", max.cutoff = "q95",
ncol = 3, reduction = "tsne_adt")FeaturePlot(pbmc, features = features.adt, min.cutoff = "q05", max.cutoff = "q95",
ncol = 3, reduction = "tsne_rna")FeaturePlot(pbmc, features = features.hto, min.cutoff = "q05", max.cutoff = "q95",
ncol = 3, reduction = "tsne_rna")